% 定义 Newton 插值函数
function P = newton_interpolation(x, y, xx)
    n = length(x);
    a = divided_difference(x, y);
    P = a(n) * ones(size(xx));
    for k = n-1:-1:1
        P = a(k) + (xx - x(k)).*P;
    end
end

% 定义除法差分的函数
function a = divided_difference(x, y)
    n = length(y);
    v = zeros(n,n);
    v(:,1) = y(:); 
    for j=2:n
        for k=j:n
            v(k,j) = (v(k,j-1)-v(k-1,j-1))/(x(k)-x(k-j+1));
        end
    end
    a = diag(v);
end
